home *** CD-ROM | disk | FTP | other *** search
/ Linux Cubed Series 3: Developer Tools / Linux Cubed Series 3 - Developer Tools.iso / devel / lang / forth / pfe-0.000 / pfe-0 / pfe-0.9.13 / src / dblsub.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-07-17  |  5.7 KB  |  259 lines
  1. /*
  2.  * This file is part of the portable Forth environment written in ANSI C.
  3.  * Copyright (C) 1995  Dirk Uwe Zoller
  4.  *
  5.  * This library is free software; you can redistribute it and/or
  6.  * modify it under the terms of the GNU Library General Public
  7.  * License as published by the Free Software Foundation; either
  8.  * version 2 of the License, or (at your option) any later version.
  9.  *
  10.  * This library is distributed in the hope that it will be useful,
  11.  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
  13.  * See the GNU Library General Public License for more details.
  14.  *
  15.  * You should have received a copy of the GNU Library General Public
  16.  * License along with this library; if not, write to the Free
  17.  * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  *
  19.  * This file is version 0.9.13 of 17-July-95
  20.  * Check for the latest version of this package via anonymous ftp at
  21.  *    roxi.rz.fht-mannheim.de:/pub/languages/forth/pfe-VERSION.tar.gz
  22.  * or    sunsite.unc.edu:/pub/languages/forth/pfe-VERSION.tar.gz
  23.  * or    ftp.cygnus.com:/pub/forth/pfe-VERSION.tar.gz
  24.  *
  25.  * Please direct any comments via internet to
  26.  *    duz@roxi.rz.fht-mannheim.de.
  27.  * Thank You.
  28.  */
  29. /*
  30.  * dblsub.c ---       Subroutines for double number (64 Bit) arithmetics.
  31.  * (duz 16Jul93)
  32.  */
  33.  
  34. #include "forth.h"
  35. #include "dblsub.h"
  36.  
  37. void
  38. dasl (dCell *a, int n)        /* left shift of *a by n positions */
  39. {
  40.   while (--n >= 0)
  41.     {
  42.       a->hi <<= 1;
  43.       a->hi += HIGHBIT (a->lo);
  44.       a->lo <<= 1;
  45.     }
  46. }
  47.  
  48. void
  49. dasr (dCell *a, int n)        /* arithm. right shift of *a by n positions */
  50. {
  51.   while (--n >= 0)
  52.     {
  53.       a->lo >>= 1;
  54.       a->lo += a->hi << (CELLBITS - 1);
  55.       a->hi >>= 1;
  56.     }
  57. }
  58. /* *INDENT-OFF* */
  59. void
  60. madd (dCell * a, Cell b)    /* add b to a */
  61. {
  62.   uCell c;            /* carry */
  63.  
  64.   D3 (*a) = c  = (uCell) D3 (*a) + W1 (b); c >>= HALFCELL;
  65.   D2 (*a) = c += (uCell) D2 (*a) + W0 (b); c >>= HALFCELL;
  66.   D1 (*a) = c += (uCell) D1 (*a);       c >>= HALFCELL;
  67.   D0 (*a) = c +  (uCell) D0 (*a);
  68. }
  69.  
  70. void
  71. dadd (dCell * a, dCell * b)    /* add b to a */
  72. {
  73.   uCell c;            /* carry */
  74.  
  75.   D3 (*a) = c  = (uCell) D3 (*a) + D3 (*b); c >>= HALFCELL;
  76.   D2 (*a) = c += (uCell) D2 (*a) + D2 (*b); c >>= HALFCELL;
  77.   D1 (*a) = c += (uCell) D1 (*a) + D1 (*b); c >>= HALFCELL;
  78.   D0 (*a) = c +  (uCell) D0 (*a) + D0 (*b);
  79. }
  80.  
  81. void
  82. dsub (dCell * a, dCell * b)    /* subtract b from a */
  83. {
  84.   Cell c;            /* carry */
  85.  
  86.   D3 (*a) = c  = (Cell) D3 (*a) - (Cell) D3 (*b); c >>= HALFCELL;
  87.   D2 (*a) = c += (Cell) D2 (*a) - (Cell) D2 (*b); c >>= HALFCELL;
  88.   D1 (*a) = c += (Cell) D1 (*a) - (Cell) D1 (*b); c >>= HALFCELL;
  89.   D0 (*a) = c +  (Cell) D0 (*a) - (Cell) D0 (*b);
  90. }
  91.  
  92. void
  93. dnegate (dCell * a)        /* negate a */
  94. {
  95.   Cell c;            /* carry */
  96.  
  97.   D3 (*a) = c = -(Cell) D3 (*a); c >>= HALFCELL;
  98.   D2 (*a) = c -= (Cell) D2 (*a); c >>= HALFCELL;
  99.   D1 (*a) = c -= (Cell) D1 (*a); c >>= HALFCELL;
  100.   D0 (*a) = c -  (Cell) D0 (*a);
  101. }
  102. /* *INDENT-ON */
  103.  
  104. int
  105. dless (dCell * a, dCell * b)    /* result: a < b */
  106. {
  107.   return a->hi != b->hi
  108.   ? a->hi < b->hi
  109.   : a->lo < b->lo;
  110. }
  111.  
  112. int
  113. duless (udCell * a, udCell * b)    /* result: a < b */
  114. {
  115.   return a->hi != b->hi ? a->hi < b->hi : a->lo < b->lo;
  116. }
  117.  
  118. udCell
  119. ummul (uCell a, uCell b)    /* unsigned multiply, mixed precision */
  120. {
  121.   udCell res;
  122.   uCell c, p;
  123.  
  124.   res.lo = (uCell) W1 (a) * W1 (b);
  125.   if (W0 (a))
  126.     {
  127.       p = (uCell) W0 (a) * W1 (b);
  128.       if (W0 (b))
  129.     {
  130.       uCell q = (uCell) W1 (a) * W0 (b);
  131.       res.hi = (uCell) W0 (a) * W0 (b);
  132.       /* *INDENT-OFF* */
  133.       D2 (res)  = c  = (uCell) D2 (res) + W1 (p) + W1 (q); c >>= HALFCELL;
  134.       D1 (res)  = c += (uCell) D1 (res) + W0 (p) + W0 (q); c >>= HALFCELL;
  135.       /* *INDENT-ON */
  136.       D0 (res) += c;
  137.     }
  138.       else
  139.     goto three;
  140.     }
  141.   else
  142.     {
  143.       if (W0 (b))
  144.     {
  145.       p = (uCell) W1 (a) * W0 (b);
  146.     three:
  147.       D2 (res) = c = (uCell) D2 (res) + W1 (p); c >>= HALFCELL;
  148.       D1 (res) = c + W0 (p);
  149.       D0 (res) = 0;
  150.     }
  151.       else
  152.     res.hi = 0;
  153.     }
  154.   return res;
  155. }
  156.  
  157. dCell
  158. mmul (Cell a, Cell b)        /* signed multiply, mixed precision */
  159. {
  160.   dCell res;
  161.   int s = 0;
  162.  
  163.   if (a < 0)
  164.     a = -a, s ^= 1;
  165.   if (b < 0)
  166.     b = -b, s ^= 1;
  167.   *(udCell *) & res = ummul (a, b);
  168.   if (s)
  169.     dnegate (&res);
  170.   return res;
  171. }
  172.  
  173. static void
  174. shift_subtract (udCell * u, uCell v)
  175. /* Divide unsigned double by single precision using shifts and subtracts.
  176.    Return quotient in u->lo, remainder in u->hi. */
  177. {
  178.   int i = CELLBITS, c = 0;
  179.   uCell q = 0, h = u->hi, l = u->lo;
  180.  
  181.   for (;;)
  182.     {
  183.       if (c || h >= v)
  184.     {
  185.       q++;
  186.       h -= v;
  187.     }
  188.       if (--i < 0)
  189.     break;
  190.       c = HIGHBIT (h);
  191.       h <<= 1;
  192.       h += HIGHBIT (l);
  193.       l <<= 1;
  194.       q <<= 1;
  195.     }
  196.   u->hi = h;
  197.   u->lo = q;
  198. }
  199.  
  200. udiv_t
  201. umdiv (udCell num, uCell denom)    /* unsigned divide procedure, mixed prec */
  202. {
  203.   udiv_t res;
  204.  
  205.   if (num.hi == 0)
  206.     {
  207.       res.quot = num.lo / denom;
  208.       res.rem = num.lo % denom;
  209.     }
  210.   else
  211.     {
  212.       shift_subtract (&num, denom);
  213.       res.quot = num.lo;
  214.       res.rem = num.hi;
  215.     }
  216.   return res;
  217. }
  218.  
  219. fdiv_t
  220. smdiv (dCell num, Cell denom)    /* symmetric divide procedure, mixed prec */
  221. {
  222.   fdiv_t res;
  223.   int sq = 0, sr = 0;
  224.  
  225.   if (num.hi < 0)
  226.     {
  227.       if (num.hi == -1 && (Cell) num.lo < 0)
  228.     goto simple;
  229.       dnegate (&num);
  230.       sq ^= 1;
  231.       sr ^= 1;
  232.     }
  233.   else if (num.hi == 0 && (Cell) num.lo > 0)
  234.     {
  235.     simple:
  236.       res.quot = (Cell) num.lo / denom;
  237.       res.rem  = (Cell) num.lo % denom;
  238.       return res;
  239.     }
  240.   if (denom < 0)
  241.     {
  242.       denom = -denom;
  243.       sq ^= 1;
  244.     }
  245.   shift_subtract ((udCell *) & num, denom);
  246.   res.quot = sq ? -num.lo : num.lo;
  247.   res.rem  = sr ? -num.hi : num.hi;
  248.   return res;
  249. }
  250.  
  251. fdiv_t
  252. fmdiv (dCell num, Cell denom)    /* floored divide procedure, mixed prec */
  253. {
  254.   fdiv_t res = smdiv (num, denom);
  255.   if (res.rem && (num.hi ^ denom) < 0)
  256.     res.quot--, res.rem += denom;
  257.   return res;
  258. }
  259.